// Copyright 2022 Jeroen van Nugteren

// Permission is hereby granted, free of charge, to any person obtaining a copy of this software
// and associated documentation files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:

// The above copyright notice and this permission notice shall be included in all copies or
// substantial portions of the Software.

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
// OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

// include header
#include "perimeter.hh"

// common headers
#include "rat/common/gmshfile.hh"
#include "rat/common/extra.hh"

// distmesh headers
#include "distfun.hh"

// code specific to Rat
namespace rat{namespace dm{

	// constructcor
	Perimeter::Perimeter(
		const arma::Mat<fltp>& Rn, 
		const arma::Mat<arma::uword> &n){
		assert(Rn.n_rows==2);
		assert(n.n_rows==2);
		Rn_ = Rn; n_ = n;
	}

	// constructcor
	Perimeter::Perimeter(const std::list<ShPerimeterPr>&perimeters){
		// allocate
		arma::field<arma::Mat<fltp> > Rn_fld(1,perimeters.size());
		arma::field<arma::Mat<arma::uword> > n_fld(1,perimeters.size());

		// walk oevr perimeters
		arma::uword cnt = 0; arma::uword node_shift = 0;
		for(auto it = perimeters.begin();it!=perimeters.end();it++){
			// get mesh
			Rn_fld(cnt) = (*it)->get_nodes();
			n_fld(cnt) = (*it)->get_elements() + node_shift;

			// increment counters
			node_shift+=Rn_fld(cnt).n_cols; cnt++; 
		}

		// combine
		Rn_ = cmn::Extra::field2mat(Rn_fld);
		n_ = cmn::Extra::field2mat(n_fld);
	}

	// factory
	ShPerimeterPr Perimeter::create(
		const arma::Mat<fltp>& Rn, 
		const arma::Mat<arma::uword> &n){
		return std::make_shared<Perimeter>(Rn,n);
	}

	// factory
	ShPerimeterPr Perimeter::create(
		const std::list<ShPerimeterPr>&perimeters){
		return std::make_shared<Perimeter>(perimeters);
	}
		
	// getters
	const arma::Mat<fltp>& Perimeter::get_nodes()const{
		return Rn_;
	}

	// setters
	const arma::Mat<arma::uword>& Perimeter::get_elements()const{
		return n_;
	}

	// intersection function
	arma::Mat<fltp> Perimeter::intersect_df(const ShDistFunPr& df){
		// no intersection points
		if(Rn_.empty())return arma::Mat<fltp>(0,2);

		// tolerance for gradient (if gradient is zero shapes are coinciding)
		const fltp tol = RAT_CONST(1e-6);

		// evaluate distance function at nodes
		const arma::Col<fltp> distance = df->calc_distance(Rn_.t());

		// get distance to boundary at start and end of each element
		const arma::Col<fltp> d1 = distance.rows(n_.row(0));
		const arma::Col<fltp> d2 = distance.rows(n_.row(1));

		// find elements that are on the boundary of the distance function
		const arma::Col<arma::uword> idx_root = arma::find((d1<0.0 && d2>=0.0 && (d2 - d1)>tol) || (d1>0.0 && d2<=0.0 && (d1 - d2)>tol));

		// allocate
		arma::Mat<fltp> Pi(idx_root.n_elem,2); arma::uword cnt = 0;

		// walk over roots
		for(arma::uword i=0;i<idx_root.n_elem;i++){
			// get index of element
			const arma::uword eid = idx_root(i);

			// get element nodes
			const arma::Col<arma::uword>::fixed<2> nroot = n_.col(eid);
			
			// element vector
			const arma::Col<fltp>::fixed<2> V = arma::diff(Rn_.cols(nroot),1,1);
		
			// relative distance
			const fltp dv = d2(eid) - d1(eid);

			// find intersection point (approximate)
			Pi.row(cnt++) = (Rn_.col(nroot(0)) + (-d1(eid)/dv)*V).t();
		}

		// return points
		return Pi;
	}

	// // refine intersection point
	// void Perimeter::refine_intersection(
	// 	arma::Mat<fltp> &Pi, 
	// 	const ShDistFunPr& df1, 
	// 	const ShDistFunPr& df2,
	// 	const fltp delta,
	// 	const fltp abstol){

	// 	// settings
	// 	const arma::uword max_iter = 1024;
	// 	const fltp damping = 0.4;

	// 	// finite difference steps
	// 	arma::Row<fltp>::fixed<2> dx(arma::fill::zeros); dx(0) = delta;
	// 	arma::Row<fltp>::fixed<2> dy(arma::fill::zeros); dy(1) = delta;

	// 	// iterations
	// 	for(arma::uword i=0;i<max_iter;i++){
	// 		// check distance against tolerance
	// 		if(arma::all(arma::abs(df1->calc_distance(Pi))<abstol && 
	// 			arma::abs(df2->calc_distance(Pi))<abstol))break;

	// 		std::cout<<df1->calc_distance(Pi)<<std::endl;
	// 		std::cout<<df2->calc_distance(Pi)<<std::endl;

	// 		// alternate
	// 		const ShDistFunPr& df = i%2==0 ? df1 : df2;
	// 		const ShDistFunPr& df_other = i%2==1 ? df1 : df2;
		
	// 		// calculate distance function
	// 		const arma::Col<fltp> d0 = df->calc_distance(Pi);

	// 		// calculate finite difference points
	// 		const arma::Col<fltp> d1 = df->calc_distance(Pi.each_row()-dx/2);
	// 		const arma::Col<fltp> d2 = df->calc_distance(Pi.each_row()+dx/2);
	// 		const arma::Col<fltp> d3 = df->calc_distance(Pi.each_row()-dy/2);
	// 		const arma::Col<fltp> d4 = df->calc_distance(Pi.each_row()+dy/2);

	// 		// finite difference
	// 		const arma::Col<fltp> gradx = (d2 - d1)/delta;
	// 		const arma::Col<fltp> grady = (d4 - d3)/delta;
	// 		const arma::Col<fltp> dgrad = arma::hypot(gradx,grady);

	// 		// step
	// 		const arma::Mat<fltp> dd = arma::join_horiz(d0%gradx/dgrad,d0%grady/dgrad);

	// 		// find points that are not gradient zero
	// 		arma::Col<arma::uword> idx_grad = arma::find(dgrad!=0);

	// 		// bring point to nearest boundary
	// 		Pi.rows(idx_grad) -= damping*dd.rows(idx_grad);

	// 		// now search along boundary
	// 		const arma::Col<fltp> a = df_other->calc_distance(Pi);
	// 		const arma::Col<fltp> b = df_other->calc_distance(Pi + delta*arma::join_horiz(-grady/dgrad, gradx/dgrad));
	// 		const arma::Col<fltp> grad2 = (b-a)/delta;
	// 		const arma::Mat<fltp> dd2 = (a/grad2)%arma::join_horiz(-grady/dgrad, gradx/dgrad).eval().each_col();
	// 		arma::Col<arma::uword> idx_grad2 = arma::find(grad2!=0 && dgrad!=0);
	// 		Pi.rows(idx_grad2) -= damping*dd2.rows(idx_grad2);

	// 		// check
	// 		if(i==max_iter-1)rat_throw_line("could not solve position of fixed point");
	// 	}
	// }

	// // refine intersection point
	// void Perimeter::drop_lowgradangle(
	// 	arma::Mat<fltp> &Pi, 
	// 	const ShDistFunPr& df1, 
	// 	const ShDistFunPr& df2,
	// 	const fltp delta){

	// 	// finite difference steps
	// 	arma::Row<fltp>::fixed<2> dx(arma::fill::zeros); dx(0) = delta;
	// 	arma::Row<fltp>::fixed<2> dy(arma::fill::zeros); dy(1) = delta;

	// 	// calculate finite difference points
	// 	const arma::Col<fltp> d11 = df1->calc_distance(Pi.each_row()-dx/2);
	// 	const arma::Col<fltp> d12 = df1->calc_distance(Pi.each_row()+dx/2);
	// 	const arma::Col<fltp> d13 = df1->calc_distance(Pi.each_row()-dy/2);
	// 	const arma::Col<fltp> d14 = df1->calc_distance(Pi.each_row()+dy/2);
	// 	const arma::Col<fltp> d21 = df2->calc_distance(Pi.each_row()-dx/2);
	// 	const arma::Col<fltp> d22 = df2->calc_distance(Pi.each_row()+dx/2);
	// 	const arma::Col<fltp> d23 = df2->calc_distance(Pi.each_row()-dy/2);
	// 	const arma::Col<fltp> d24 = df2->calc_distance(Pi.each_row()+dy/2);

	// 	// finite difference
	// 	const arma::Col<fltp> gradx1 = (d12 - d11)/delta;
	// 	const arma::Col<fltp> grady1 = (d14 - d13)/delta;
	// 	const arma::Col<fltp> gradx2 = (d22 - d21)/delta;
	// 	const arma::Col<fltp> grady2 = (d24 - d23)/delta;

	// 	// calculate norm
	// 	const arma::Col<fltp> dgrad1 = arma::hypot(gradx1,grady1);
	// 	const arma::Col<fltp> dgrad2 = arma::hypot(gradx2,grady2);

	// 	// calculate angle
	// 	const arma::Col<fltp> angle = arma::acos((gradx1%gradx2 + grady1%grady2)/(dgrad1%dgrad2));

	// 	// drop rows that have small angle between gradients
	// 	Pi.shed_rows(arma::find(angle<30*arma::Datum<fltp>::tau/360));
	// }



	

}}
